home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1993 July / InfoMagic USENET CD-ROM July 1993.ISO / sources / misc / volume11 / starchart / part18 < prev    next >
Encoding:
Text File  |  1990-03-25  |  32.2 KB  |  1,110 lines
  1. Newsgroups: comp.sources.misc
  2. subject: v11i046: starchart 3.2 Part 18/32
  3. from: ccount@ATHENA.MIT.EDU
  4. Sender: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc)
  5.  
  6. Posting-number: Volume 11, Issue 46
  7. Submitted-by: ccount@ATHENA.MIT.EDU
  8. Archive-name: starchart/part18
  9.  
  10. #! /bin/sh
  11. # This is a shell archive.  Remove anything before this line, then unpack
  12. # it by saving it into a file and typing "sh file".  To overwrite existing
  13. # files, type "sh file -c".  You can also feed this as standard input via
  14. # unshar, or by typing "sh <file", e.g..  If this archive is complete, you
  15. # will see the following message at the end:
  16. #        "End of archive 18 (of 32)."
  17. # Contents:  starchart/ssup.c.ab
  18. PATH=/bin:/usr/bin:/usr/ucb ; export PATH
  19. if test -f 'starchart/ssup.c.ab' -a "${1}" != "-c" ; then 
  20.   echo shar: Will not clobber existing file \"'starchart/ssup.c.ab'\"
  21. else
  22. echo shar: Extracting \"'starchart/ssup.c.ab'\" \(30231 characters\)
  23. sed "s/^X//" >'starchart/ssup.c.ab' <<'END_OF_FILE'
  24. X  double xroot1, xroot2;
  25. X  double yr1, yr2, r1, r2;
  26. X  int n;
  27. X  int x1, y1;
  28. X  double lat_1, lon_1;
  29. X  int in;
  30. X
  31. X  if (fabs(x_2 - x_1) < 1e-5) {        /* Line has infinite slope */
  32. X    *x = x_1;
  33. X    *y = a*x_1*x_1 + b;
  34. X    *int_1 = TRUE;
  35. X  } else {                /* Line slope may be calculated */
  36. X    c = (y_2 - y_1)/(x_2 - x_1);
  37. X    d = y_1 - c * x_1;
  38. X    if (a < 1e-5) {            /* virtually a straight line y = b */
  39. X      if (fabs(c) < 1e-5) {        /* Constant y */
  40. X    n = 0;
  41. X      } else {
  42. X    xroot1 = (b - d)/c;
  43. X    n = 1;
  44. X      };
  45. X    } else {
  46. X      quadrat(a, -c, b-d, &xroot1, &xroot2, &n);
  47. X    };
  48. X    if (n == 0) {            /* No intersection */
  49. X      *int_1 = FALSE;
  50. X    } else if (n == 1) {        /* One intersection */
  51. X      *x = xroot1;
  52. X      *y = a*xroot1*xroot1 + b;
  53. X      *int_1 = TRUE;
  54. X    } else {                /* Two intersections */
  55. X      yr1 = c*xroot1 + d;
  56. X      yr2 = c*xroot2 + d;
  57. X      r1 = (xroot1-x_1)*(xroot1-x_1) + (yr1-y_1)*(yr1-y_1)
  58. X        + (xroot1-x_2)*(xroot1-x_2) + (yr1-y_2)*(yr1-y_2);
  59. X      r2 = (xroot2-x_1)*(xroot2-x_1) + (yr2-y_1)*(yr2-y_1)
  60. X        + (xroot2-x_2)*(xroot2-x_2) + (yr2-y_2)*(yr2-y_2);
  61. X      if (r1 > r2) {
  62. X    *x = xroot2;
  63. X    *y = yr2;
  64. X    *int_1 = TRUE;
  65. X      } else {
  66. X    *x = xroot1;
  67. X    *y = yr1;
  68. X    *int_1 = TRUE;
  69. X      };
  70. X    }
  71. X  }
  72. X
  73. X  if (*int_1)
  74. X    if (((((y_1-Fuz) <= *y) && (*y <= (y_2+Fuz)))
  75. X     || (((y_2-Fuz) <= *y) && (*y <= (y_1+Fuz))))
  76. X    && ((((x_1-Fuz) <= *x) && (*x <= (x_2+Fuz)))
  77. X        || (((x_2-Fuz) <= *x) && (*x <= (x_1+Fuz)))))
  78. X      *int_1 = TRUE;
  79. X    else
  80. X      *int_1 = FALSE;
  81. X
  82. X  if (*int_1) {
  83. X    inv_gt(*x, *y, &lat_1, &lon_1);
  84. X    xform(lat_1, lon_1, &x1, &y1, &in);
  85. X    if (!in) *int_1 = FALSE;
  86. X  }
  87. X}
  88. X
  89. X
  90. Xellip_intersect(x_1, y_1, x_2, y_2, a, b, y0, x, y, int_1)
  91. X     double x_1, y_1, x_2, y_2, *x, *y;
  92. X     double a, b, y0;        /* x*x/a + (y-y0)*(y-y0)/b - 1 = 0 */
  93. X     int *int_1;
  94. X{
  95. X  double c, d;
  96. X  double xroot1, xroot2;
  97. X  double yr1, yr2, r1, r2;
  98. X  int n;
  99. X  int x1, y1;
  100. X  double lat_1, lon_1;
  101. X  int in;
  102. X
  103. X  if (fabs(x_2 - x_1) < 1e-5) {        /* Line has infinite slope */
  104. X    xroot1 = xroot2 = *x = x_1;
  105. X    if (x_1*x_1/a > 1.0) {
  106. X      n = 0;
  107. X      *int_1 = FALSE;
  108. X    } else if (x_1*x_1/a == 1.0) {
  109. X      yr1 = y0;
  110. X      n = 1;
  111. X      *int_1 = TRUE;
  112. X    } else {
  113. X      yr1 = y0 + sqrt(b*(1 - x_1*x_1/a));
  114. X      yr2 = y0 - sqrt(b*(1 - x_1*x_1/a));
  115. X      n = 2;
  116. X      *int_1 = TRUE;
  117. X    }
  118. X  } else {                /* Line slope may be calculated */
  119. X    c = (y_2 - y_1)/(x_2 - x_1);
  120. X    d = y_1 - c * x_1;
  121. X    quadrat(b + a*c*c, 2*a*c*(d - y0), a*(d - y0)*(d - y0) - a*b,
  122. X        &xroot1, &xroot2, &n);
  123. X    if (n == 0) {        /* No intersection */
  124. X      *int_1 = FALSE;
  125. X    } else if (n == 1) {    /* One intersection */
  126. X      *x = xroot1;
  127. X      *y = c*xroot1 + d;
  128. X      *int_1 = TRUE;
  129. X    } else {            /* Two intersections */
  130. X      yr1 = c*xroot1 + d;
  131. X      yr2 = c*xroot2 + d;
  132. X      *int_1 = TRUE;
  133. X    };
  134. X  };
  135. X
  136. X  
  137. X  if (n == 2) {
  138. X    r1 = (xroot1-x_1)*(xroot1-x_1) + (yr1-y_1)*(yr1-y_1)
  139. X          + (xroot1-x_2)*(xroot1-x_2) + (yr1-y_2)*(yr1-y_2);
  140. X    r2 = (xroot2-x_1)*(xroot2-x_1) + (yr2-y_1)*(yr2-y_1)
  141. X          + (xroot2-x_2)*(xroot2-x_2) + (yr2-y_2)*(yr2-y_2);
  142. X    if (r1 > r2) {
  143. X      *x = xroot2;
  144. X      *y = yr2;
  145. X      *int_1 = TRUE;
  146. X    } else {
  147. X      *x = xroot1;
  148. X      *y = yr1;
  149. X      *int_1 = TRUE;
  150. X    }
  151. X  }
  152. X
  153. X  if (*int_1)
  154. X    if ((((y_1 <= *y) && (*y <= y_2)) || ((y_2 <= *y) && (*y <= y_1)))
  155. X    && (((x_1 <= *x) && (*x <= x_2)) || ((x_2 <= *x) && (*x <= x_1))))
  156. X      *int_1 = TRUE;
  157. X    else
  158. X      *int_1 = FALSE;
  159. X
  160. X  if (*int_1) {
  161. X    inv_gt(*x, *y, &lat_1, &lon_1);
  162. X    xform(lat_1, lon_1, &x1, &y1, &in);
  163. X    if (!in) *int_1 = FALSE;
  164. X  }
  165. X}
  166. X
  167. X
  168. Xcirc_intersect(x_1, y_1, x_2, y_2, r, x1, y1, int_1, x2, y2, int_2)
  169. X     double x_1, y_1, x_2, y_2, *x1, *y1, *x2, *y2;
  170. X     double r;
  171. X     int *int_1, *int_2;
  172. X{
  173. X  double c, d;
  174. X  double xroot1, xroot2;
  175. X  double yr1, yr2, r1, r2;
  176. X  int n;
  177. X  int xt1, yt1;
  178. X  double lat_1, lon_1;
  179. X  int in;
  180. X
  181. X  if (fabs(x_2 - x_1) < 1e-5) {        /* Line has infinite slope */
  182. X    xroot1 = xroot2 = *x1 = *x2 = x_1;
  183. X    if (fabs(r) > fabs(x_1)) {
  184. X      yr1 = sqrt(r*r - x_1*x_1);
  185. X      yr2 = -yr1;
  186. X      n = 2;
  187. X      *int_1 = *int_2 = TRUE;
  188. X    } else if (fabs(r) == fabs(x_1)) {
  189. X      yr1 = 0;
  190. X      n = 1;
  191. X      *int_1 = *int_2 = TRUE;
  192. X    } else {
  193. X      n = 0;
  194. X      *int_1 = *int_2 = FALSE;
  195. X    };
  196. X  } else {                /* Line slope may be calculated */
  197. X    c = (y_2 - y_1)/(x_2 - x_1);
  198. X    d = y_1 - c * x_1;
  199. X    quadrat(1 + c*c, 2*c*d, d*d - r*r,
  200. X        &xroot1, &xroot2, &n);
  201. X    if (n == 0) {        /* No intersection */
  202. X      *int_1 = *int_2 = FALSE;
  203. X    } else if (n == 1) {    /* One intersection */
  204. X      *x1 = xroot1;
  205. X      *y1 = c*xroot1 + d;
  206. X      *int_1 = TRUE;
  207. X      *int_2 = FALSE;
  208. X    } else {            /* Two intersections */
  209. X      yr1 = c*xroot1 + d;
  210. X      yr2 = c*xroot2 + d;
  211. X      *int_1 = *int_2 = TRUE;
  212. X    };
  213. X  };
  214. X
  215. X  if (n == 2) {
  216. X    r1 = (xroot1-x_1)*(xroot1-x_1) + (yr1-y_1)*(yr1-y_1)
  217. X          + (xroot1-x_2)*(xroot1-x_2) + (yr1-y_2)*(yr1-y_2);
  218. X    r2 = (xroot2-x_1)*(xroot2-x_1) + (yr2-y_1)*(yr2-y_1)
  219. X          + (xroot2-x_2)*(xroot2-x_2) + (yr2-y_2)*(yr2-y_2);
  220. X    if (r1 > r2) {
  221. X      *x1 = xroot2;
  222. X      *y1 = yr2;
  223. X      *x2 = xroot1;
  224. X      *y2 = yr1;
  225. X      *int_1 = *int_2 = TRUE;
  226. X    } else {
  227. X      *x1 = xroot1;
  228. X      *y1 = yr1;
  229. X      *x2 = xroot2;
  230. X      *y2 = yr2;
  231. X      *int_1 = *int_2 = TRUE;
  232. X    }
  233. X  }
  234. X
  235. X  if (*int_1)
  236. X    if ((((y_1 <= *y1) && (*y1 <= y_2)) || ((y_2 <= *y1) && (*y1 <= y_1)))
  237. X    && (((x_1 <= *x1) && (*x1 <= x_2)) || ((x_2 <= *x1) && (*x1 <= x_1))))
  238. X      *int_1 = TRUE;
  239. X    else
  240. X      *int_1 = FALSE;
  241. X
  242. X  if (*int_1) {
  243. X    inv_gt(*x1, *y1, &lat_1, &lon_1);
  244. X    xform(lat_1, lon_1, &xt1, &yt1, &in);
  245. X    if (!in) *int_1 = FALSE;
  246. X  }
  247. X
  248. X  if (*int_2)
  249. X    if ((((y_1 <= *y2) && (*y2 <= y_2)) || ((y_2 <= *y2) && (*y2 <= y_1)))
  250. X    && (((x_1 <= *x2) && (*x2 <= x_2)) || ((x_2 <= *x2) && (*x2 <= x_1))))
  251. X      *int_2 = TRUE;
  252. X    else
  253. X      *int_2 = FALSE;
  254. X
  255. X  if (*int_2) {
  256. X    inv_gt(*x2, *y2, &lat_1, &lon_1);
  257. X    xform(lat_1, lon_1, &xt1, &yt1, &in);
  258. X    if (!in) *int_2 = FALSE;
  259. X  }
  260. X}
  261. X
  262. X
  263. Xquadrat(a, b, c, x_1, x_2, n)
  264. X     double a, b, c, *x_1, *x_2;
  265. X     int *n;
  266. X{
  267. X  double t;
  268. X
  269. X  if (a == 0) {
  270. X    *n = 0;
  271. X  } else {
  272. X    t = b * b - 4 * a * c;
  273. X    if (t < 0) {
  274. X      *n = 0;
  275. X    } else if (t == 0) {
  276. X      *x_1 = -b/(2*a);
  277. X      *n = 1;
  278. X    } else {
  279. X      *x_1 = (-b + sqrt(t))/(2*a);
  280. X      *x_2 = (-b - sqrt(t))/(2*a);
  281. X      *n = 2;
  282. X    };
  283. X  };
  284. X}
  285. X
  286. X/* Draw a curved line between points 1 and 2.
  287. Xclipxform has been called, xloc1, yloc1, xloc2, yloc2 are in bounds */
  288. X
  289. Xdrawcurveline(lat1, lon1, lat2, lon2,
  290. X          xloc1, yloc1, xloc2, yloc2, line_style, great_circle, clevel)
  291. X     int xloc1, yloc1, xloc2, yloc2;
  292. X     double lat1, lon1, lat2, lon2;
  293. X     int line_style;
  294. X     int great_circle;     /* draw as great circle */
  295. X     int clevel;         /* safety valve */
  296. X{
  297. X  double mlat, mlon;    /* midpoint lat and long */
  298. X  double tlat1, tlat2;    /* temporary */
  299. X  int txloc1, tyloc1, in;
  300. X  int txloc2, tyloc2;
  301. X  double slope1;
  302. X  int mxloc, myloc;    /* transformed */
  303. X  int mpx, mpy;        /* from given x,y */
  304. X  int inregion;
  305. X
  306. X/* ra difference should be less than 180 degrees: take shortest path */
  307. X  if ((xfs_proj_mode == STEREOGR) || (xfs_proj_mode == GNOMONIC)
  308. X      || (xfs_proj_mode == ORTHOGR))
  309. X    if ((lon1 - lon2) > 180.0) lon1 -= 360.0;
  310. X    else if ((lon2 - lon1) > 180.0) lon2 -= 360.0;
  311. X
  312. X/* Adjust so lon1 and lon2 are continuous across region: see xform */
  313. X/* needed so midpoint is correct */
  314. X  if ((xfs_proj_mode == SANSONS) || (xfs_proj_mode == RECTANGULAR)) {
  315. X    if ((xf_west < 0.0) && (lon1>(xf_west+360.0))) lon1 -= 360.0;
  316. X    if ((xf_east > 360.0) && (lon1<(xf_east-360.0))) lon1 += 360.0;
  317. X    if ((xf_west < 0.0) && (lon2>(xf_west+360.0))) lon2 -= 360.0;
  318. X    if ((xf_east > 360.0) && (lon2<(xf_east-360.0))) lon2 += 360.0;
  319. X
  320. X    /* path crosses boundary of map */
  321. X    /* if the distance from point1 to left (East) + point2 to right
  322. X       or point2 to right and point2 to left
  323. X       is less than distance from point1 to point2,
  324. X       then we have a problem. */
  325. X    if (xfs_wide_warn)
  326. X      if ((fabs(lon1-xf_east) + fabs(xf_west-lon2)) < (fabs(lon1-lon2))) {
  327. X    slope1 = (lat2 - lat1)/(fabs(lon1-xf_east) + fabs(xf_west-lon2));
  328. X
  329. X    tlat1 = lat1 + slope1 * fabs(lon1-(xf_east-xf_c_scale));
  330. X    xform(tlat1, xf_east-xf_c_scale, &txloc1, &tyloc1, &in);
  331. X    drawcurveline(lat1, lon1, tlat1, xf_east-xf_c_scale, xloc1, yloc1,
  332. X              txloc1, tyloc1, line_style, great_circle, 0);
  333. X
  334. X    tlat2 = lat2 - slope1 * fabs(lon2-(xf_west+xf_c_scale));
  335. X    xform(tlat2, xf_west+xf_c_scale, &txloc2, &tyloc2, &in);
  336. X    drawcurveline(tlat2, xf_west+xf_c_scale, lat2, lon2, txloc2, tyloc2,
  337. X              xloc2, yloc2, line_style, great_circle, 0);
  338. X    return;
  339. X      } else if ((fabs(lon2-xf_east)+fabs(xf_west-lon1)) < (fabs(lon1-lon2))) {
  340. X    slope1 = (lat1 - lat2)/(fabs(lon2-xf_east) + fabs(xf_west-lon1));
  341. X
  342. X    tlat1 = lat2 + slope1 * fabs(lon2-(xf_east-xf_c_scale));
  343. X    xform(tlat1, xf_east-xf_c_scale, &txloc1, &tyloc1, &in);
  344. X    drawcurveline(lat2, lon2, tlat1, xf_east-xf_c_scale, xloc2, yloc2, 
  345. X              txloc1, tyloc1, line_style, great_circle, 0);
  346. X
  347. X    tlat2 = lat1 - slope1 * fabs(lon1-(xf_west+xf_c_scale));
  348. X    xform(tlat2, xf_west+xf_c_scale, &txloc2, &tyloc2, &in);
  349. X    drawcurveline(tlat2, xf_west+xf_c_scale, lat1, lon1, txloc2, tyloc2,
  350. X              xloc1, yloc1, line_style, great_circle, 0);
  351. X    return;
  352. X      };
  353. X  }
  354. X
  355. X
  356. X  if (great_circle) {
  357. X    gcmidpoint(lat1, lon1, lat2, lon2, &mlat, &mlon);
  358. X  } else {
  359. X    mlat = (lat1 + lat2)/2;
  360. X    mlon = (lon1 + lon2)/2;
  361. X  };
  362. X
  363. X
  364. X
  365. X  xform(mlat, mlon, &mxloc, &myloc, &inregion);
  366. X
  367. X  mpx = (xloc1 + xloc2)/2;
  368. X  mpy = (yloc1 + yloc2)/2;
  369. X
  370. X  if (((abs(mpx - mxloc) + abs(mpy - myloc)) > 0) && (clevel < 100)) {
  371. X    /* center is not where it should be */
  372. X    drawcurveline(lat1, lon1, mlat, mlon, xloc1, yloc1,
  373. X          mxloc, myloc, line_style, great_circle, ++clevel);
  374. X    drawcurveline(mlat, mlon, lat2, lon2, mxloc, myloc,
  375. X          xloc2, yloc2, line_style, great_circle, ++clevel);
  376. X  } else {
  377. X    D_movedraw(xloc1, yloc1, xloc2, yloc2, line_style);
  378. X  }
  379. X}
  380. X
  381. Xgcmidpoint(lat1, lon1, lat2, lon2, pmlat, pmlon)
  382. X     double lat1, lon1, lat2, lon2, *pmlat, *pmlon;
  383. X{
  384. X  double l1, m1, n1;
  385. X  double l2, m2, n2;
  386. X  double l3, m3, n3;
  387. X
  388. X  /* transform from (ra, dec) to l m n direction cosines */
  389. X  l1 = DCOS(lat1)*DCOS(lon1);
  390. X  m1 = DCOS(lat1)*DSIN(lon1);
  391. X  n1 = DSIN(lat1);
  392. X
  393. X  l2 = DCOS(lat2)*DCOS(lon2);
  394. X  m2 = DCOS(lat2)*DSIN(lon2);
  395. X  n2 = DSIN(lat2);
  396. X
  397. X  l3 = l1 + l2;
  398. X  m3 = m1 + m2;
  399. X  n3 = n1 + n2;
  400. X  n3 /= sqrt(l3*l3 + m3*m3 + n3*n3);
  401. X
  402. X  *pmlon = atan2(m3, l3) / 0.0174532925199;
  403. X  if ((*pmlon < 0) && (lon1 > 0) && (lon2 > 0)) *pmlon += 360.0;
  404. X  *pmlat = asin(n3) / 0.0174532925199;
  405. X}
  406. X
  407. X
  408. X
  409. X/* Area handling */
  410. X/* Works often, but is far from rigorous */
  411. Xstatic struct t_area_struct {
  412. X  double lat, lon;
  413. X  int great_circle;
  414. X  int clipped_at;
  415. X} area_raw[100], area_clip[200], area_1[200], area_2[200];
  416. Xstatic int nsegments = 0, nseg_clip = 0, nseg_1 = 0, nseg_2 = 0;
  417. X
  418. Xareastart(lat, lon, great_circle)
  419. X     double lat, lon;
  420. X     int great_circle;
  421. X{
  422. X  if (nsegments != 0)
  423. X    areafinish();
  424. X
  425. X  nsegments = 0;
  426. X  area_raw[nsegments].lat = lat;
  427. X  area_raw[nsegments].lon = lon;
  428. X  area_raw[nsegments].great_circle = great_circle;
  429. X  nsegments++;
  430. X}
  431. X
  432. X
  433. Xareaadd(lat, lon, great_circle)
  434. X     double lat, lon;
  435. X     int great_circle;
  436. X{
  437. X  area_raw[nsegments].lat = lat;
  438. X  area_raw[nsegments].lon = lon;
  439. X  area_raw[nsegments].great_circle = great_circle;
  440. X  nsegments++;
  441. X}
  442. X
  443. Xint wrap_area();
  444. X
  445. Xareafinish()
  446. X{
  447. X  int i;
  448. X  int xloc, yloc, xloc2, yloc2;
  449. X  int inregion;
  450. X  double tlat1, tlon1, tlat2, tlon2, tlat3, tlon3, tlat4, tlon4;
  451. X  double slope1;
  452. X  int done, old_n;
  453. X  int curr_area;
  454. X  int gt12;
  455. X  int wraps;
  456. X  int done_clip;
  457. X  done_clip = FALSE;
  458. X  nseg_clip = nseg_1 = nseg_2 = 0;
  459. X  wraps = FALSE;
  460. X
  461. X
  462. X  /* skip if no point is within region */
  463. X  done = TRUE;
  464. X  for (i = 0; i < nsegments; i++) {
  465. X    xform(area_raw[i].lat, area_raw[i].lon, &xloc, &yloc, &inregion);
  466. X    if (inregion) done = FALSE;
  467. X  };
  468. X  if (done) return;
  469. X
  470. X  if ((xfs_proj_mode == SANSONS) || (xfs_proj_mode == RECTANGULAR)) {
  471. X    if (xfs_wide_warn)
  472. X      for (i = 0; i < (nsegments - 1); i++)
  473. X    if (wrap_area(area_raw[i].lat, area_raw[i].lon,
  474. X              area_raw[i+1].lat, area_raw[i+1].lon, 
  475. X              &tlat1, &tlon1, &tlat2, &tlon2))
  476. X      wraps = TRUE;
  477. X
  478. X    if (!xfs_wide_warn) {
  479. X      /* Use algorithm to clip area_raw to rectangular region */
  480. X      /* Use area_1 and area_2 as scratch */
  481. X
  482. X      /* clip against west from area_raw to area_1 */
  483. X      tlat1 = area_raw[nsegments-1].lat;
  484. X      tlon1 = area_raw[nsegments-1].lon;
  485. X      for (i = 0, nseg_1 = 0; i < nsegments; i++) {
  486. X    tlat4 = area_raw[i].lat;
  487. X    tlon4 = area_raw[i].lon;
  488. X
  489. X    if (eastof(tlon4, xf_west)) {
  490. X      if (eastof(tlon1, xf_west)) {
  491. X        /* last point and new point (4) are in, add new point */
  492. X      } else {
  493. X        /* new point in, but old out.
  494. X         Calculate intersection point (2), add it and new point */
  495. X        slope1 = (tlat4 - tlat1) / (tlon4 - tlon1);
  496. X
  497. X        tlat2 = tlat1 + slope1 * (xf_west-tlon1);
  498. X        tlon2 = xf_west;
  499. X        area_1[nseg_1].lat = tlat2;
  500. X        area_1[nseg_1].lon = tlon2;
  501. X        area_1[nseg_1].great_circle = area_raw[i].great_circle;
  502. X        area_1[nseg_1].clipped_at = WEST_CLIP;
  503. X        nseg_1++;
  504. X      }
  505. X      area_1[nseg_1].lat = tlat4;
  506. X      area_1[nseg_1].lon = tlon4;
  507. X      area_1[nseg_1].great_circle = area_raw[i].great_circle;
  508. X      area_1[nseg_1].clipped_at = 0;
  509. X      nseg_1++;
  510. X    } else {
  511. X      if (eastof(tlon1, xf_west)) {
  512. X        /* new point out, but old point in.
  513. X           Calculate intersection and add it */
  514. X        slope1 = (tlat4 - tlat1) / (tlon4 - tlon1);
  515. X
  516. X        tlat2 = tlat1 + slope1 * (xf_west-tlon1);
  517. X        tlon2 = xf_west;
  518. X        area_1[nseg_1].lat = tlat2;
  519. X        area_1[nseg_1].lon = tlon2;
  520. X        area_1[nseg_1].great_circle = area_raw[i].great_circle;
  521. X        area_1[nseg_1].clipped_at = WEST_CLIP;
  522. X        nseg_1++;
  523. X      }
  524. X    }
  525. X    tlat1 = tlat4;
  526. X    tlon1 = tlon4;
  527. X      }
  528. X
  529. X
  530. X      /* clip against south from area_1 to area_2 */
  531. X      tlat1 = area_1[nseg_1-1].lat;
  532. X      tlon1 = area_1[nseg_1-1].lon;
  533. X      for (i = 0, nseg_2 = 0; i < nseg_1; i++) {
  534. X    tlat4 = area_1[i].lat;
  535. X    tlon4 = area_1[i].lon;
  536. X
  537. X    if (tlat4 > xf_south) {
  538. X      if (tlat1 > xf_south) {
  539. X        /* last point and new point (4) are in, add new point */
  540. X      } else {
  541. X        /* new point in, but old out.
  542. X         Calculate intersection point (2), add it and new point */
  543. X        slope1 = (tlon4 - tlon1) / (tlat4 - tlat1);
  544. X
  545. X        tlat2 = xf_south;
  546. X        tlon2 = tlon1 + slope1 * (xf_south-tlat1);
  547. X        area_2[nseg_2].lat = tlat2;
  548. X        area_2[nseg_2].lon = tlon2;
  549. X        area_2[nseg_2].great_circle = area_1[i].great_circle;
  550. X        area_2[nseg_2].clipped_at = SOUTH_CLIP;
  551. X        nseg_2++;
  552. X      }
  553. X      area_2[nseg_2].lat = tlat4;
  554. X      area_2[nseg_2].lon = tlon4;
  555. X      area_2[nseg_2].great_circle = area_1[i].great_circle;
  556. X      area_2[nseg_2].clipped_at = 0;
  557. X      nseg_2++;
  558. X    } else {
  559. X      if (tlat1 > xf_south) {
  560. X        /* new point out, but old point in.
  561. X           Calculate intersection and add it */
  562. X        slope1 = (tlon4 - tlon1) / (tlat4 - tlat1);
  563. X
  564. X        tlat2 = xf_south;
  565. X        tlon2 = tlon1 + slope1 * (xf_south-tlat1);
  566. X        area_2[nseg_2].lat = tlat2;
  567. X        area_2[nseg_2].lon = tlon2;
  568. X        area_2[nseg_2].great_circle = area_1[i].great_circle;
  569. X        area_2[nseg_2].clipped_at = SOUTH_CLIP;
  570. X        nseg_2++;
  571. X      }
  572. X    }
  573. X    tlat1 = tlat4;
  574. X    tlon1 = tlon4;
  575. X      }
  576. X
  577. X      /* clip against east from area_2 to area_1 */
  578. X      tlat1 = area_2[nseg_2-1].lat;
  579. X      tlon1 = area_2[nseg_2-1].lon;
  580. X      for (i = 0, nseg_1 = 0; i < nseg_2; i++) {
  581. X    tlat4 = area_2[i].lat;
  582. X    tlon4 = area_2[i].lon;
  583. X
  584. X    if (westof(tlon4, xf_east)) {
  585. X      if (westof(tlon1, xf_east)) {
  586. X        /* last point and new point (4) are in, add new point */
  587. X      } else {
  588. X        /* new point in, but old out.
  589. X         Calculate intersection point (2), add it and new point */
  590. X        slope1 = (tlat4 - tlat1) / (tlon4 - tlon1);
  591. X
  592. X        tlat2 = tlat1 + slope1 * (xf_east-tlon1);
  593. X        tlon2 = xf_east;
  594. X        area_1[nseg_1].lat = tlat2;
  595. X        area_1[nseg_1].lon = tlon2;
  596. X        area_1[nseg_1].great_circle = area_2[i].great_circle;
  597. X        area_1[nseg_1].clipped_at = EAST_CLIP;
  598. X        nseg_1++;
  599. X      }
  600. X      area_1[nseg_1].lat = tlat4;
  601. X      area_1[nseg_1].lon = tlon4;
  602. X      area_1[nseg_1].great_circle = area_2[i].great_circle;
  603. X      area_1[nseg_1].clipped_at = 0;
  604. X      nseg_1++;
  605. X    } else {
  606. X      if (westof(tlon1, xf_east)) {
  607. X        /* new point out, but old point in.
  608. X           Calculate intersection and add it */
  609. X        slope1 = (tlat4 - tlat1) / (tlon4 - tlon1);
  610. X
  611. X        tlat2 = tlat1 + slope1 * (xf_east-tlon1);
  612. X        tlon2 = xf_east;
  613. X        area_1[nseg_1].lat = tlat2;
  614. X        area_1[nseg_1].lon = tlon2;
  615. X        area_1[nseg_1].great_circle = area_2[i].great_circle;
  616. X        area_1[nseg_1].clipped_at = EAST_CLIP;
  617. X        nseg_1++;
  618. X      }
  619. X    }
  620. X    tlat1 = tlat4;
  621. X    tlon1 = tlon4;
  622. X      }
  623. X
  624. X      /* clip against north from area_1 to area_clip */
  625. X      tlat1 = area_1[nseg_1-1].lat;
  626. X      tlon1 = area_1[nseg_1-1].lon;
  627. X      for (i = 0, nseg_clip = 0; i < nseg_1; i++) {
  628. X    tlat4 = area_1[i].lat;
  629. X    tlon4 = area_1[i].lon;
  630. X
  631. X    if (tlat4 < xf_north) {
  632. X      if (tlat1 < xf_north) {
  633. X        /* last point and new point (4) are in, add new point */
  634. X      } else {
  635. X        /* new point in, but old out.
  636. X         Calculate intersection point (2), add it and new point */
  637. X        slope1 = (tlon4 - tlon1) / (tlat4 - tlat1);
  638. X
  639. X        tlat2 = xf_north;
  640. X        tlon2 = tlon1 + slope1 * (xf_north-tlat1);
  641. X        area_clip[nseg_clip].lat = tlat2;
  642. X        area_clip[nseg_clip].lon = tlon2;
  643. X        area_clip[nseg_clip].great_circle = area_1[i].great_circle;
  644. X        area_clip[nseg_clip].clipped_at = NORTH_CLIP;
  645. X        nseg_clip++;
  646. X      }
  647. X      area_clip[nseg_clip].lat = tlat4;
  648. X      area_clip[nseg_clip].lon = tlon4;
  649. X      area_clip[nseg_clip].great_circle = area_1[i].great_circle;
  650. X      area_clip[nseg_clip].clipped_at = 0;
  651. X      nseg_clip++;
  652. X    } else {
  653. X      if (tlat1 < xf_north) {
  654. X        /* new point out, but old point in.
  655. X           Calculate intersection and add it */
  656. X        slope1 = (tlon4 - tlon1) / (tlat4 - tlat1);
  657. X
  658. X        tlat2 = xf_north;
  659. X        tlon2 = tlon1 + slope1 * (xf_north-tlat1);
  660. X        area_clip[nseg_clip].lat = tlat2;
  661. X        area_clip[nseg_clip].lon = tlon2;
  662. X        area_clip[nseg_clip].great_circle = area_1[i].great_circle;
  663. X        area_clip[nseg_clip].clipped_at = NORTH_CLIP;
  664. X        nseg_clip++;
  665. X      }
  666. X    }
  667. X    tlat1 = tlat4;
  668. X    tlon1 = tlon4;
  669. X      }
  670. X
  671. X      if ((area_clip[0].lat != area_clip[nseg_clip-1].lat)
  672. X      || (area_clip[0].lon != area_clip[nseg_clip-1].lon)) {
  673. X    area_clip[nseg_clip].lat = area_clip[0].lat;
  674. X    area_clip[nseg_clip].lon = area_clip[0].lon;
  675. X    area_clip[nseg_clip].great_circle = area_raw[i].great_circle;
  676. X    area_clip[nseg_clip].clipped_at = WEST_CLIP;
  677. X    nseg_clip++;
  678. X      }
  679. X
  680. X      if (nseg_clip > 0) doarea(area_clip, nseg_clip);
  681. X
  682. X      done_clip = TRUE;
  683. X    };
  684. X  };
  685. X
  686. X  if (!done_clip) {
  687. X    /* Scan list of segments, constructing clipped polygon */
  688. X    for (i = 0; i < (nsegments - 1); i++) {
  689. X      if (clipr_xform(area_raw[i].lat, area_raw[i].lon,
  690. X              area_raw[i+1].lat, area_raw[i+1].lon, 
  691. X              &xloc, &yloc, &xloc2, &yloc2, area_raw[i+1].great_circle,
  692. X              &tlat1, &tlon1, &tlat2, &tlon2)) {
  693. X    addsegment(tlat1, tlon1, tlat2, tlon2,
  694. X           (nseg_clip > 0) ? area_clip[nseg_clip-1].clipped_at : 0,
  695. X           clip_at1, clip_at2,
  696. X           area_raw[i].great_circle, area_raw[i+1].great_circle);
  697. X      };
  698. X    };
  699. X
  700. X    /* ensure that area is closed. */
  701. X    /* redraw first segment */
  702. X    if ((nseg_clip > 0) && ((area_clip[nseg_clip-1].lat != area_clip[0].lat)
  703. X               || (area_clip[nseg_clip-1].lon != area_clip[0].lon))) {
  704. X      i = 0;
  705. X      done = FALSE;
  706. X      old_n = nseg_clip;
  707. X      do {
  708. X    if (clipr_xform(area_raw[i].lat, area_raw[i].lon,
  709. X            area_raw[i+1].lat, area_raw[i+1].lon, 
  710. X            &xloc, &yloc, &xloc2, &yloc2,
  711. X            area_raw[i+1].great_circle,
  712. X            &tlat1, &tlon1, &tlat2, &tlon2)) {
  713. X      addsegment(tlat1, tlon1, tlat2, tlon2,
  714. X             (nseg_clip > 0) ? area_clip[nseg_clip-1].clipped_at : 0,
  715. X             clip_at1, clip_at2,
  716. X             area_raw[i].great_circle, area_raw[i+1].great_circle);
  717. X      done = TRUE;
  718. X    }
  719. X    i++;
  720. X      } while ((i < (nsegments - 1)) && (!done));
  721. X      if ((old_n+1) < nseg_clip) nseg_clip--; /* only needed one added */
  722. X    }
  723. X
  724. X
  725. X    /* Now fix problem of areas crossing from left to right across boundary */
  726. X    /* if the distance from point1 to left (East) + point2 to right
  727. X       or point2 to right and point2 to left
  728. X       is less than distance from point1 to point2,
  729. X       then we have a problem. */
  730. X    curr_area = 1;
  731. X    for (i = 0; i < (nseg_clip - 1); i++) {
  732. X      if (wraps && wrap_area(area_clip[i].lat, area_clip[i].lon,
  733. X                 area_clip[i+1].lat, area_clip[i+1].lon, 
  734. X                 &tlat1, &tlon1, &tlat2, &tlon2)) {
  735. X    /* Crosses boundary */
  736. X    switch (curr_area) {
  737. X    case 1:
  738. X                /* finish segment in part 1,
  739. X                   begin part 2 */
  740. X      area_1[nseg_1].lat = area_clip[i].lat;
  741. X      area_1[nseg_1].lon = area_clip[i].lon;
  742. X      area_1[nseg_1].great_circle = area_clip[i].great_circle;
  743. X      area_1[nseg_1].clipped_at = area_clip[i].clipped_at;
  744. X      nseg_1++;
  745. X
  746. X      area_1[nseg_1].lat = tlat1;
  747. X      area_1[nseg_1].lon = tlon1;
  748. X      area_1[nseg_1].great_circle = area_clip[i].great_circle;
  749. X      area_1[nseg_1].clipped_at = area_clip[i].clipped_at;
  750. X      nseg_1++;
  751. X
  752. X      curr_area = 2;
  753. X
  754. X      area_2[nseg_2].lat = tlat2;
  755. X      area_2[nseg_2].lon = tlon2;
  756. X      area_2[nseg_2].great_circle = area_clip[i].great_circle;
  757. X      area_2[nseg_2].clipped_at = area_clip[i].clipped_at;
  758. X      nseg_2++;
  759. X
  760. X      gt12 = (tlon1 > tlon2); /* area 1 greater than area 2 */
  761. X      break;
  762. X    case 2:
  763. X                /* finish part 2, draw it,
  764. X                   close gap in part 1 */
  765. X      if (((gt12) && (tlon1 < tlon2))
  766. X          || ((!gt12) && (tlon1 > tlon2))){
  767. X        tlat3 = tlat1;
  768. X        tlon3 = tlon1;
  769. X        tlat1 = tlat2;
  770. X        tlon1 = tlon2;
  771. X        tlat2 = tlat3;
  772. X        tlon2 = tlon3;
  773. X      };
  774. X      area_2[nseg_2].lat = area_clip[i].lat;
  775. X      area_2[nseg_2].lon = area_clip[i].lon;
  776. X      area_2[nseg_2].great_circle = area_clip[i].great_circle;
  777. X      area_2[nseg_2].clipped_at = area_clip[i].clipped_at;
  778. X      nseg_2++;
  779. X
  780. X      area_2[nseg_2].lat = tlat2;
  781. X      area_2[nseg_2].lon = tlon2;
  782. X      area_2[nseg_2].great_circle = area_clip[i].great_circle;
  783. X      area_2[nseg_2].clipped_at = area_clip[i].clipped_at;
  784. X      nseg_2++;
  785. X
  786. X      area_2[nseg_2].lat = area_2[0].lat;
  787. X      area_2[nseg_2].lon = area_2[0].lon;
  788. X      area_2[nseg_2].great_circle = TRUE;
  789. X      area_2[nseg_2].clipped_at = area_2[0].clipped_at;
  790. X      nseg_2++;
  791. X
  792. X      doarea(area_2, nseg_2);
  793. X
  794. X      curr_area = 1;
  795. X
  796. X      area_1[nseg_1].lat = tlat1;
  797. X      area_1[nseg_1].lon = tlon1;
  798. X      area_1[nseg_1].great_circle = TRUE;
  799. X      area_1[nseg_1].clipped_at = area_1[nseg_1-1].clipped_at;
  800. X      nseg_1++;
  801. X
  802. X      break;
  803. X    default:
  804. X      break;
  805. X    };
  806. X      } else {
  807. X    switch (curr_area) {
  808. X    case 1:
  809. X      area_1[nseg_1].lat = area_clip[i].lat;
  810. X      area_1[nseg_1].lon = area_clip[i].lon;
  811. X      area_1[nseg_1].great_circle = area_clip[i].great_circle;
  812. X      area_1[nseg_1].clipped_at = area_clip[i].clipped_at;
  813. X      nseg_1++;
  814. X      break;
  815. X    case 2:
  816. X      area_2[nseg_2].lat = area_clip[i].lat;
  817. X      area_2[nseg_2].lon = area_clip[i].lon;
  818. X      area_2[nseg_2].great_circle = area_clip[i].great_circle;
  819. X      area_2[nseg_2].clipped_at = area_clip[i].clipped_at;
  820. X      nseg_2++;
  821. X      break;
  822. X    default:
  823. X      break;
  824. X    }
  825. X      }
  826. X    };
  827. X    if (nseg_clip > 0) {
  828. X      area_1[nseg_1].lat = area_clip[nseg_clip-1].lat;
  829. X      area_1[nseg_1].lon = area_clip[nseg_clip-1].lon;
  830. X      area_1[nseg_1].great_circle = area_clip[nseg_clip-1].great_circle;
  831. X      area_1[nseg_1].clipped_at = area_clip[nseg_clip-1].clipped_at;
  832. X      nseg_1++;
  833. X    };
  834. X
  835. X    if (nseg_1 > 0) doarea(area_1, nseg_1);
  836. X  };
  837. X
  838. X  /* Mark it done */
  839. X  nsegments = 0;
  840. X}
  841. X
  842. Xdoarea(area, nseg)
  843. X     struct t_area_struct area[];
  844. X     int nseg;
  845. X{
  846. X  int i;
  847. X  int xloc, yloc, xloc2, yloc2;
  848. X  double tlat1, tlon1, tlat2, tlon2;
  849. X
  850. X  for (i = 0; i < (nseg - 1); i++) {
  851. X    if (clipr_xform(area[i].lat, area[i].lon,
  852. X            area[i+1].lat, area[i+1].lon, 
  853. X            &xloc, &yloc, &xloc2, &yloc2, area[i+1].great_circle,
  854. X            &tlat1, &tlon1, &tlat2, &tlon2)) {
  855. X
  856. X      if (i == 0) D_areamove(xloc, yloc);
  857. X      /* Here we can add tests to add segments along the boundary of
  858. X     projection modes other than SANSONS and RECTANGULAR */
  859. X
  860. X      addareacurved(tlat1, tlon1, tlat2, tlon2,
  861. X            xloc, yloc, xloc2, yloc2,
  862. X            area[i+1].great_circle, 0);
  863. X      if (i == (nseg - 2)) D_areafill(xloc2, yloc2);
  864. X    }
  865. X  }
  866. X}
  867. X
  868. Xint wrap_area(lat1, lon1, lat2, lon2, plat1, plon1, plat2, plon2)
  869. X     double lat1, lon1, lat2, lon2;
  870. X     double *plat1, *plon1, *plat2, *plon2;
  871. X{
  872. X  double slope1, tlat1, tlat2;
  873. X
  874. X/* Adjust so lon1 and lon2 are continuous across region: see xform */
  875. X/* needed so midpoint is correct */
  876. X  if ((xfs_proj_mode == SANSONS) || (xfs_proj_mode == RECTANGULAR)) {
  877. X    if ((xf_west < 0.0) && (lon1>(xf_west+360.0))) lon1 -= 360.0;
  878. X    if ((xf_east > 360.0) && (lon1<(xf_east-360.0))) lon1 += 360.0;
  879. X    if ((xf_west < 0.0) && (lon2>(xf_west+360.0))) lon2 -= 360.0;
  880. X    if ((xf_east > 360.0) && (lon2<(xf_east-360.0))) lon2 += 360.0;
  881. X
  882. X    /* path crosses boundary of map */
  883. X    /* if the distance from point1 to left (East) + point2 to right
  884. X       or point2 to right and point2 to left
  885. X       is less than distance from point1 to point2,
  886. X       then we have a problem. */
  887. X    if (xfs_wide_warn)
  888. X      if ((fabs(lon1-xf_east) + fabs(xf_west-lon2)) < (fabs(lon1-lon2))) {
  889. X    /* point 1 close to east boundary */
  890. X    slope1 = (lat2 - lat1)/(fabs(lon1-xf_east) + fabs(xf_west-lon2));
  891. X
  892. X    tlat1 = lat1 + slope1 * fabs(lon1-(xf_east-xf_c_scale));
  893. X
  894. X    /* segment from (lat1, lon1) to (tlat1, xf_east-xf_c_scale)
  895. X       is on one side of boundary */
  896. X    *plat1 = tlat1;
  897. X    *plon1 = xf_east-xf_c_scale;
  898. X    if (*plon1 > 360) *plon1 -= 360;
  899. X    if (*plon1 < 0) *plon1 += 360;
  900. X
  901. X    tlat2 = lat2 - slope1 * fabs(lon2-(xf_west+xf_c_scale));
  902. X    /* segment from (tlat2, xf_west+xf_c_scale) to (lat2, lon2)
  903. X       is on one side of boundary */
  904. X    *plat2 = tlat2;
  905. X    *plon2 = xf_west+xf_c_scale;
  906. X    if (*plon2 > 360) *plon2 -= 360;
  907. X    if (*plon2 < 0) *plon2 += 360;
  908. X    return TRUE;
  909. X      } else if ((fabs(lon2-xf_east)+fabs(xf_west-lon1)) < (fabs(lon1-lon2))) {
  910. X    /* point 1 close to west boundary */
  911. X    slope1 = (lat1 - lat2)/(fabs(lon2-xf_east) + fabs(xf_west-lon1));
  912. X
  913. X    tlat1 = lat2 + slope1 * fabs(lon2-(xf_east-xf_c_scale));
  914. X
  915. X    /* segment from (lat2, lon2) to (tlat1, xf_east-xf_c_scale)
  916. X       is on one side of boundary */
  917. X    *plat2 = tlat1;
  918. X    *plon2 = xf_east-xf_c_scale;
  919. X    if (*plon2 > 360) *plon2 -= 360;
  920. X    if (*plon2 < 0) *plon2 += 360;
  921. X
  922. X    tlat2 = lat1 - slope1 * fabs(lon1-(xf_west+xf_c_scale));
  923. X    /* segment from (tlat2, xf_west+xf_c_scale) to (lat1, lon1)
  924. X       is on one side of boundary */
  925. X    *plat1 = tlat2;
  926. X    *plon1 = xf_west+xf_c_scale;
  927. X    if (*plon1 > 360) *plon1 -= 360;
  928. X    if (*plon1 < 0) *plon1 += 360;
  929. X    return TRUE;
  930. X      };
  931. X  }
  932. X  return FALSE;
  933. X}
  934. X
  935. X
  936. Xaddsegment(tlat1, tlon1, tlat2, tlon2, last_clip, clip_1, clip_2,
  937. X       gc_1, gc_2)
  938. X     double tlat1, tlon1, tlat2, tlon2;
  939. X     int last_clip, clip_1, clip_2;
  940. X     int gc_1, gc_2;
  941. X{
  942. X
  943. X  /* just add a segment */
  944. X  area_clip[nseg_clip].lat = tlat1;
  945. X  area_clip[nseg_clip].lon = tlon1;
  946. X  area_clip[nseg_clip].great_circle = gc_1;
  947. X  area_clip[nseg_clip].clipped_at = clip_1;
  948. X  nseg_clip++;
  949. X  area_clip[nseg_clip].lat = tlat2;
  950. X  area_clip[nseg_clip].lon = tlon2;
  951. X  area_clip[nseg_clip].great_circle = gc_2;
  952. X  area_clip[nseg_clip].clipped_at = clip_2;
  953. X  nseg_clip++;
  954. X}
  955. X
  956. X
  957. X
  958. X/* addareacurved: add a set of boundary segments between points 1 and 2
  959. X   clipxform has been called, xloc1, yloc1, xloc2, yloc2 are in bounds */
  960. Xaddareacurved(lat1, lon1, lat2, lon2,
  961. X          xloc1, yloc1, xloc2, yloc2, great_circle, clevel)
  962. X     int xloc1, yloc1, xloc2, yloc2;
  963. X     double lat1, lon1, lat2, lon2;
  964. X     int great_circle;     /* draw as great circle */
  965. X     int clevel;         /* safety valve */
  966. X{
  967. X  double mlat, mlon;    /* midpoint lat and long */
  968. X  int mxloc, myloc;    /* transformed */
  969. X  int mpx, mpy;        /* from given x,y */
  970. X  int inregion;
  971. X  
  972. X/* ra difference should be less than 180 degrees: take shortest path */
  973. X  if ((xfs_proj_mode == STEREOGR) || (xfs_proj_mode == GNOMONIC)
  974. X      || (xfs_proj_mode == ORTHOGR))
  975. X    if ((lon1 - lon2) > 180.0) lon1 -= 360.0;
  976. X    else if ((lon2 - lon1) > 180.0) lon2 -= 360.0;
  977. X
  978. X/* Adjust so lon1 and lon2 are continuous across region: see xform */
  979. X/* needed so midpoint is correct */
  980. X  if ((xfs_proj_mode == SANSONS) || (xfs_proj_mode == RECTANGULAR)) {
  981. X    if ((xf_west < 0.0) && (lon1>(xf_west+360.0))) lon1 -= 360.0;
  982. X    if ((xf_east > 360.0) && (lon1<(xf_east-360.0))) lon1 += 360.0;
  983. X    if ((xf_west < 0.0) && (lon2>(xf_west+360.0))) lon2 -= 360.0;
  984. X    if ((xf_east > 360.0) && (lon2<(xf_east-360.0))) lon2 += 360.0;
  985. X  }
  986. X
  987. X
  988. X  if (great_circle) {
  989. X    gcmidpoint(lat1, lon1, lat2, lon2, &mlat, &mlon);
  990. X  } else {
  991. X    mlat = (lat1 + lat2)/2;
  992. X    mlon = (lon1 + lon2)/2;
  993. X  };
  994. X
  995. X  xform(mlat, mlon, &mxloc, &myloc, &inregion);
  996. X
  997. X  mpx = (xloc1 + xloc2)/2;
  998. X  mpy = (yloc1 + yloc2)/2;
  999. X
  1000. X  if (((abs(mpx - mxloc) + abs(mpy - myloc)) > 0) && (clevel < 100)) {
  1001. X    /* center is not where it should be */
  1002. X    addareacurved(lat1, lon1, mlat, mlon,
  1003. X          xloc1, yloc1, mxloc, myloc, great_circle, ++clevel);
  1004. X    addareacurved(mlat, mlon, lat2, lon2,
  1005. X          mxloc, myloc, xloc2, yloc2, great_circle, ++clevel);
  1006. X  } else {
  1007. X    D_areaadd(xloc2, yloc2);
  1008. X  }
  1009. X}
  1010. X
  1011. X
  1012. X
  1013. X
  1014. X/* Translate two digit encoded size string */
  1015. X/* Maximum parsed: 89000 secs of arc = 24.666 degrees */
  1016. X/* Warning, should use 32 bit integers to allow this value */
  1017. Xlong size_obj(size_str)
  1018. X     char *size_str;
  1019. X{
  1020. X  char chr1, chr2;
  1021. X
  1022. X  chr1 = islower(size_str[0]) ? toupper(size_str[0]) : size_str[0];
  1023. X  chr2 = islower(size_str[1]) ? toupper(size_str[1]) : size_str[1];
  1024. X
  1025. X  if ((chr1 == ' ') && (chr2 == ' ')) return -1;
  1026. X  if (chr1 == ' ') return chr2 - '0';
  1027. X  if (isdigit(chr1)) return ((chr1-'0')*10L + (chr2-'0'));
  1028. X  if (chr1 < 'J') return ((chr1-'A'+1)*100L + (chr2-'0')*10L);
  1029. X  if (chr1 < 'S') return ((chr1-'I')*1000L + (chr2-'0')*100L);
  1030. X  if (chr1 <= 'Z') return ((chr1-'R')*10000L + (chr2-'0')*1000L);
  1031. X  return -1;
  1032. X}
  1033. X
  1034. X/*
  1035. XExamples:
  1036. X
  1037. X"  "        -1
  1038. X" 6"        6
  1039. X"09"        9
  1040. X"73"        73
  1041. X"A0"        100
  1042. X"C3"        330
  1043. X"D5"        450
  1044. X"I6"        960
  1045. X"J2"        1200
  1046. X"R3"        9300
  1047. X"S6"        16000
  1048. X"Z0"        80000
  1049. X"Z9"        89000
  1050. X*/
  1051. X
  1052. X/* Precession functions */
  1053. X/* precession based on formulae from
  1054. X   Astronomical Almanac, 1988, p B19
  1055. X*/
  1056. Xstatic double M, N;
  1057. Xinitprecess(eout)
  1058. X     double eout;
  1059. X{
  1060. X  double T;
  1061. X
  1062. X  T = (eout - 2000.0) / 100.0;
  1063. X  M = (1.2812323 + (0.0003879 + 0.0000101 * T) * T) * T;
  1064. X  N = (0.5567530 - (0.0001185 - 0.0000116 * T) * T) * T;
  1065. X}
  1066. X
  1067. Xdo_precess(lat, lon)
  1068. X     double *lat, *lon;
  1069. X{
  1070. X  double am, dm, a2, d2;
  1071. X
  1072. X  /* am = /alpha_m, dm = /delta_m */
  1073. X
  1074. X  am = *lon + (M + N * DSIN(*lon) * DTAN(*lat))/2.0;
  1075. X  dm = *lat + N*DCOS(am)/2.0;
  1076. X  a2 = *lon + M + N*DSIN(am)*DTAN(dm);
  1077. X  d2 = *lat + N * DCOS(am);
  1078. X
  1079. X  if (a2 >= 360.0) a2 -= 360.0;
  1080. X  if (a2 < 0.0) a2 += 360.0;
  1081. X
  1082. X  *lon = a2;
  1083. X  *lat = d2;
  1084. X}
  1085. X
  1086. END_OF_FILE
  1087. if test 30231 -ne `wc -c <'starchart/ssup.c.ab'`; then
  1088.     echo shar: \"'starchart/ssup.c.ab'\" unpacked with wrong size!
  1089. fi
  1090. # end of 'starchart/ssup.c.ab'
  1091. fi
  1092. echo shar: End of archive 18 \(of 32\).
  1093. cp /dev/null ark18isdone
  1094. MISSING=""
  1095. for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 ; do
  1096.     if test ! -f ark${I}isdone ; then
  1097.     MISSING="${MISSING} ${I}"
  1098.     fi
  1099. done
  1100. if test "${MISSING}" = "" ; then
  1101.     echo You have unpacked all 32 archives.
  1102.     rm -f ark[1-9]isdone ark[1-9][0-9]isdone
  1103. else
  1104.     echo You still need to unpack the following archives:
  1105.     echo "        " ${MISSING}
  1106. fi
  1107. ##  End of shell archive.
  1108. exit 0
  1109.  
  1110.